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ABSTRACT 

We use deep adaptive mesh refinement simulations of isothermal self-gravitating supersonic turbulence to 
study the imprints of gravity on the mass density distribution in molecular clouds. The simulations show that 
the density distribution in self-gravitating clouds develops an extended power-law tail at high densities on top 
of the usual lognormal. We associate the origin of the tail with self-similar collapse solutions and predict 
the power index values in the range from -7/4 to -3/2 that agree with both simulations and observations of 
star- forming molecular clouds. 
Subject headings: ISM: structure — methods: numerical — stars: formation — turbulence 



1. INTRODUCTION 

The probability density function (PDF) of the mass density 
in non-self-gravitating isother mal superso nic turbulence 
is believed to be lognorm a l dVazguez-Se madeni 1994; 
Passot & Vazquez-Semadenil 119981; iPadoan & Nordlundl 
1999; Ki-itsuk et al. 2007^ Some hints of power-law 
tails, however, have been noticed in numerical simula- 
tions of the self-gravitating turbulent interstellar medium 



(ISM ; e.g.. IScal o et al. 1998; Klessen 2000; Dib & Burkert 
2005'; 'Slvz et al. 2005; Vazauez-SemadenieLaLl i2008t 
Feden-ath et al. 2008; Collins et al. 2010). More recently 
similar tails were also found in high dyna mic range obser- 
vations of star-formin g mol ecular clouds (.Kainulainen et al.l 
120091; lLombardiet"anl2010h . While it is understood that the 
density PDF holds the key to phenomenology of star forma- 
tion (e^g^Pado an & Nordlund 2002; Ki-umholz & McKee 
I2OOI iHermebelle & Chabried 120081: |Cho & Kim ;2011l) . the 
origin of the power-law tail in self-gravitating supersonic 
turbulence still awaits a credible explanation. A related 
question pertains to the power index value for the tail. 

Slvz et al. (2005) find a slope of -1.5 in non-m agnetic kpc- 
scale i nterstellar turbulence simulations, while iCollins et al.l 
(I2OIOI) measured -1 .6 in isothermal adaptive mesh refinement 
(AMR) MHD simulations of supersonic molecular cloud tur- 
bulence in a 10 pc box. Is there a universal power index value 
that applies to self-gravitating isothermal turbulence? What 
determines the slope? To address these questions, we ana- 
lyze a deep AMR simulation with a linear dynamic range of 
5x10^ that follows the star formation process from turbulent 
initial conditions on a scale of a few pc down to a few AU. 

We describe the simulation detail in the following section, 
while Section 3 presents the analysis of the density distribu- 
tion and provides testable predictions for the power index val- 
ues. Section 4 discusses the limitations of the model and ef- 
fects of the magnetic field on the density PDF. The final sec- 
tion outlines conclusions. 

2. NUMERICAL EXPERIMENT 

Our star formation simulation was perfor med with the 
ENZO code for cosmology and astrophysics (O'She aetaLl 
|2005) and discussed earlier in Padoan et al. (2005). We 
solve the hydrodynamic equations, including self-gravity and 
a large-scale random force to drive the turbulence. We also 
adopt an isothermal equation of state and periodic boundary 



conditions. In this simulation, AMR is automatically carried 
out in collapsing regions in order to properly resolve the Jeans 
length (Truelove et al. 1997). We use root grid of 512^ cells 
and five AMR levels with a refinement factor of four. The 
computational box has a size L = 5 pc, and the gravitational 
collapse of dense protostellar cores is resolved down to the 
scale of 2 AU. The temperature is uniform, T = 10 K, and the 
sound speed is constant, Cs = 0.2 km s~'. The mean density 
«o(H2) = 500 cm"-' and the rms flow velocity of 1.1km s"', 
typical of molecular clouds on scales of ~ 5 pc, correspond to 
sonic Mach number M^ sa 6. The free-fall time 
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the dynamical time 

Myn 



' 37r 
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1.6 Myr, 



2.3 Myr, 



and the virial parameter for this model 
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(2) 



(3) 



correspond to a 3.4 x 10^ Mq molecular cloud prone to col- 
lapse on its free-fall timescale. Indeed, the dendrogram anal- 
ysis apphed to a s napshot from a larger 1024^ simulation 
jKritsuk et a l. 2007) resembling the initial conditions adopted 
here, indicated the presence of gravitationally bound objects 
on essentially all scales within the computational domain 
(iRosolowsky et al.l l2008'). 

We began the simulation as a uniform grid turbulence 
model by stirring the gas in the computational domain for 
4.8rdyn with a large-scale random force that includes 40% di- 
latational and 60% solenoidal power and then, at f = 0, turned 
the forcing off to continue the simulation with AMR and self- 
gravity for about 0.29rdyn ~ 0.43%. 

3. EFFECTS OF SELF-GRAVITY 

Figure [T| shows the density distributions in this simula- 
tion. The red line corresponds to the initial condition at 
f = 0, when we turn on self-gravity. The density PDF at this 
time can be perfectl y represented by a lognormal distribution 
(iKritsuk et al.ll2007h . Once gravity starts to operate, a power- 
law tail develops at the high end of the distribution. After 



KRITSUK, NORMAN & WAGNER 



o 

D- 



-10 



-12 



-14 




-16 

2 4 

logiop/po 

Fig. 1. — Probability distribution functions for tlie mass density from 
an AMR simulation of self-gravitating isothermal turbulence. The red line 
shows the initial conditions corresponding to a driven, statistically steady 
Mach 6 turbulence with no self-gravity. The green Une shows the PDF for 
self-gravitating gas after 0.26% of evolution from the initial conditions. The 
effective linear dynamic range of the simulation at this instance in time is 
2048. The blue line shows a time-average PDF at 0.42 ±0.01%; time averag- 
ing helps to reduce the statistical noise at high densities. The dynamic range 
is 5 X lO'. The initial conditions can be approximated by a lognormal distri- 
bution (dashed line). A power law with a slope of —1 .695 ± 0.002 (solid line) 
provides the best fit to the high-density tail at p/po G [10, 10']. The break in 
the power index at p/po f^ lO' marks a transition to rotationally supported 
cores (slope -1, dotted line). 



0.26fff, when the creation of first-level AMR subgrids is trig- 
gered by the first collapsing objects, the density distribution 
no longer remains lognormal. As the collapse of these first 
objects proceeds, followed by further grid refinement, an ex- 
tended power-law tail emerges with a slope of about -1 .7. The 
tail departs from the initial lognormal distribution already at 
p/ Pi) ^10 and continues straight for nearly 10 dex in proba- 
bility and more than 6 dex in density. As the simulation pro- 
gresses, the slope continues to evolve slowly toward shallower 
values. The power index at the end of the simulation is -1 .67. 
An even shallower tail at very high densities, p/ po > 10^, de- 
velops as an indication of mass pile-up due to an additional 
support against gravity that comes from the conservation of 
angular momentum. ' The power index for this centrifugally 
supported part of the density distribution is very close to -1 .0. 
The fact that the power law breaks at a density slightly in ex- 
cess of 10^ po may indicate the minimum grid resolution re- 
quirement for convergence in star formation simulations with 
sink particles (roughly 32,000 for this set of parameters), but 
the main focus of this Letter is on the origin of the extended 
power law at densities below lO^po. Why does it cover over 
six orders of magnitude in mass density without a tiny bit of 
slope change? What fundamental physics is involved? 

Let us first recall that supersonic turbulence is a multi- 
scale phenomenon shaping the structure of th e ma s s dis- 
tribution in molecular cl o uds (e.g.. von Weizsacked 119511: 
iBiglari & Diamondl Il988t iKritsuk et all l2006h . At Mach 
numbers, Ms > 3, the turbulence creates a "fractal" den- 
sity distribution with the mass dimension of D^ ^2.3 
dElmegreen & Falga rond 119961: iKritsuk et all 1200 7). Since 
this mass dimension is larger than the critical value for grav- 
itational instability, D^ > D^n = 2, these highly inhomo- 

' It may partly be also due to our enforced limit on the number of al- 
lowed refinement levels (five levels maximum), which eventually violates the 



geneous system s are still subject to gravitational collapse 
( lPerdangll990h . 

The extent of the power law we obtain in the deep AMR 
simulation hints at the tail origin in the hierarchical nature of 
gravitational collapse of dense structures in molecular clouds. 
Let us take a look at the very bottom of the hierarchy, where 
AMR resolves collapsed protostellar cores. We used density 
fields for approximately cubic subvolumes centered on sev- 
eral selected dense cores to obtain the PDFs in the immediate 
vicinity of these objects. The linear size of these subvolumes 
is of order 0.005 pc. Figure |2] offers three example PDFs and 
volumetric rendering of the corresponding cores. Stretches 
of the power-law distributions are clearly present in all three 
cases, although the slopes vary from as shallow as -1.25 to 
as steep as -1.75.^ When the contributions from individual 
cores combine to form the density PDF for the whole compu- 
tational domain, by some magic the resulting slope appears to 
be the same as that from the collapsing larger-scale structures 
characterized by a lower density, p/po E [10^, 10^]. 

The easiest way to solve this puzzle is to assume that a self- 
similar collapse solution would act as a strong attractor de- 
termining the form of the density PDF in hierarchical, turbu- 
lent, self-gravitating molecular clouds. There is a large inven- 
tory of (semi-)analytical solutions for the collapse of spher- 
ically symmetric isothermal configurations to choose from. 
Whitworth & Summers (1985) arranged these similarity so- 
lutions into a banded two- dimensional continuum emb racing 
the limiting cases o f fast (JLarsonlll969t lPenstonlll969l here- 
after L P) and slow (|Shi]|ll977h collapse and their generaliza- 
tion by Hunter (191'f). While this family of solutions de- 
scribes gravitational condensation starting from a diverse set 
of idealized initial conditions, they have two important fea- 
tures in common. First, during the early stages of the evolu- 
tion preceding the formation of a singularity at the center, all 
of them develop a. p ^ r~^ density profile.^ Second, at late 
stages an expansion wave (EW) forms and propagates from 
the central singularity through the accreting material leaving 
behind a self-similar distribution with p ^ f^l'^ at r — )• 0. 

It can be readily shown that the mass density PDF for 
a spherically symmetric configuration with a p = po('"/'"o)~" 
density profile is a power law 



dV = -T^r^^d 



-■i/n 



ocdip"') 



(4) 



ITruelove et alj 119971) numerical stability condition. 



with an index m = -3/n. Thus, formally, the similarity solu- 
tions generate power-law PDFs with a slope otlp = -3/2 cor- 
responding to the r~^ profile during the early collapse stages 
and with a combination of slopes ;«lp = -3/2 (r"^ profile at 
lower densities) and wew = -2 (r"^/^ profile at higher den- 
sities) after the singularity has formed at the center. If the 
spherical symmetry is broken, for instance due to the presence 
of rotation, the situation becomes substantially more involved 
and hardly tractable analytically. Whitworth et al. (1996) sug- 
gest that a weak inward propagating compression wave that 

^ Note that the rotation-induced pile-ups at highest densities are only visi- 
ble in the two distributions that continue beyond p/po = 10'", while the third 
case shows only a hint of the pile-up at p/po > 10'. The first two cores have 
already developed relatively thin centrifugally supported disks, as can be seen 
in the renderings that show both face-on and edge-on views of the disks. The 
third core displays a rather modest flattening in the edge-on projection, indi- 
cating a weak rotation. 

' The IShiJ 09771) singular isothermal sphere (SIS) solution represents a 



hydrostatic p ' 



■ configuration from the outset. 
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Fig. 2. — Volumetric rendering (top thi'ee panels) and PDFs (bottom) of the mass density for a sample of dense cores from the AMR simulation ait = 0.42%. 
The data represent ~ 350' extractions at the highest grid resolution achieved in the simulation (2 AU). The maximum densities are pmax/po = 8.5 X lO'", 
3.2 X lO"', and 6.0 x 10** for the cores from left to right, respectively. 



triggers the formation of the singularity would converge in- 
coherently on the center and interfere with any reflected EW. 
The interference would then cause a significant degradation 
in the density profile at small radii. We believe this low- 
ers expectations for the pure EW solution in realistic situa- 
tions. However, the r~^ density profile and the correspond- 
ing slope of the density PDF otlp = -3/2 could potentially be 
preserved if an equilibrium singular disk with a flat rotation 
curve would form at the center (Norman et al. 1980; Toomre 
[l982; Havashi et al. 1982). Nevertheless, the disks formed in 
our simulation (Figure |2]i typically show steeper density pro- 
files (p ^ r~^) and rotation curves that peak at ^ 10 AU and 
then monotonically decline with radius up to 7? ^ 100 AU. 
Higher resolution would be required to follow the formation 
and structure of these disks with AMR properly. 

The PDF slope m = -1 .67 at the end of the simulation con- 
tinues to evolve slowly toward shallower values, but may still 
remain steeper than -1 .5 corresponding to the r~^ density pro- 
file. Neither does it approach -2 of the EW solution. Interest- 
ingly, a simil a r slop e of -1 .64 was independently obtained by 
IColUns et all (1201 Oh in a driven super- Alfvenic AMR MHD 
turbulence simulation with Ms = 9. This might indicate that 
some physics is missing in our discussion above. Indeed, both 
simulations model self-gravitating supersonic turbulent flows 
capable of creating the initial conditions for very dynamic 
collapses involving masses much in excess of the Bonnor- 
Ebert critical mass (lBonnodll956l: lEbertill955h . Such situa- 
tions will be better approximated by pressure-free (PF) col- 
lapse solutions (Shu 1977). The final stages of self-similar 
spherically symmetric PF collapse prior to the formation of 
central sing ularity are char acterized by the p ^ r"'^/^ den- 
sity profile (lPenstonlll969h which corresponds to the PDF 
power-law tail with nipp = -7/4. The free-fall collapse ap- 
proximation, however, always breaks down near the center 
where the effects of pressure inevitably become important, so 
the slope of the high end of the PDF should still converge to 
'Wlp = -3/2. We indeed observe a slope change from -1.7 
to -1.5 at p/po w 10^'^, see Figure [T] Since the ratio of 
gravitational-to-pressure forces in the isothermal PF collapse 
J* c>c r^/^ (lPenstonlll969l) . it scales with the mass density as 



^-1/6 -pjjjg \yeak dependence of 7* on the density is consistent 
with the appearance of break in slope at very high densities. 

The projected density of an infinite sphere with the p ^ r~" 
density distribution. 



Ys{R) = 2 



i(V^ 



^ + x^\dx(xR 



\-n 



also has a power-law PDF, 



dS ccd 



(--) 



cx£/(E''), 



(5) 



(6) 



but with a slope p = -2/(n- 1). For the LP, PF, and EW simi- 
larity solutions, p = -2, -2.8, and -4, respectively. 

Figure [3] shows the projected mass density PDFs from the 
AMR simulation at f = (red line) and t = 0.43% (blue line). 
The blue line is based on a subset of the AMR data up to 
8 AU in resolution; only the first four levels of fine mesh were 
used to make the plot. Similar to the three-dimensional den- 
sity PDF, the initial distribution is well represented by a log- 
normal. The evolved self-gravitating configuration shows a 
clear power-law tail at high column densities with a slope of 
-2.50 ±0.03. The actual slope uncertainty may be larger than 
the formally determined value of ±0.03. The distributions 
in Figure [3] show the average PDFs for all three projections, 
while only one projection would be available in real obser- 
vations. The root grid-based high-density tails for individual 
projections show some bumps (ver y similar to those seen in 
Figure 4 of Kainulai nen et al.l (2009) and some straight sec- 
tions with slopes broadly varying from as flat as -2 to steeper 
than -3. Overall, power index p = -2.5 measured in the sim- 
ulation is right in between the predicted slopes for the PF and 
LP solutions, but the EW option {p = -4) seems to be rejected. 
We expect shallower slopes in simulations with lower Mach 
numbers or in longer turbulence decay simulations without 
continuous resupply of kinetic energy from random forcing. 
Finally, a hint of a shallower slope at E > lO^*^ may indicate 
the effect of centrifugal support, similar to the slope flattening 
in Figure [T] 

4. DISCUSSION 
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Fig. 3. — Distributions of tlie projected gas density from tlie AMR simu- 
lation at ; = (red line) and at ^ = 043% (blue line). The initial distribution 
has a lognormal shape (dashed line). The final distribution has an extended 
power-law tail with a slope of —2.50 ±0.03 (solid line). 

We do not expect significant differences between driven 
and decaying turbulence models within the fra ction of the 
first free-fall time in this AMR run (see also Offner et aP 
I2OO8). By the end of the simulation, the system would still 
retain > 70% of the kinetic energy delivered by the stirring 
force, if self-gravity were not included. However, self-gravity 
slight ly increases th e kinetic energy of a turbulent system 
(e.g., Slvz et al.n 2005), thus the lack of forcing does not re- 
ally make a big difference. Since the free-fall time is shorter 
than the dynamical time that determines the energy decay 
timescale, the density PDF established by f = 0.43% can be 
only weakly sensitive to the lack of forcing. 

Since the virial parameter for this model is rather small, the 
role of self-gravity can be somewhat exaggerated. In situa- 
tions with awl one should expect a weaker effect. With a 
close to unity, assuming the same Mach number and domain 
size, the power-law tail will be shallower (more closely re- 
sembling the LP case) in a more extended density range and 
will depart from the initial lognormal distribution at densi- 
ties somewhat higher than p/po = 10. While our adopted low 
value of the virial parameter creates a distribution similar to 
that of the Taurus molecular cloud, a higher a value would 
perhaps produce a distribution more similar to that of the Lu- 
pus I cloud (see Figure 2 in Kainulainen et al. 2009). Over- 
all, it seems that cloud-to-cloud variations in the virial pa- 
rameter a and in the age of the cloud in combination with 
projection effects can account for the full diversity of high- 
density tails in the observed star-forming clouds (Figure 4 of 
iKainulainen et al. 2009). 

In this discussion on the density PDF, we so far ignored 
the effects of magnetic fields that are important for star for- 
mation. Our recent isothermal and multiphase MHD tur- 
bulence calculations, however, both show that variations in 



the level of magnetization of interstellar clouds make little 
or no difference for the h igh-density end of the density PDF 
(e.g.. iKritsuk et al]|20 10!). In the absence of self-gravity, the 
high end of the distribution preserves its lognormal shape. In 
super-Alfvenic turbulence, the magnetic field strength also 
shows a weak correlation with the gas density of the form 
B - pi/2 (iPadoan & Nordlundll999HComns et aUlMTol) . The 
most recent Zeeman splitting data fo r molecular cores also in- 
dicate a slope of 0.65 ± 0.05 (ICrutch er et al. 2010). A similar 
relation can be inferred theoretically fo r dynamically c oUaps- 
ing magnetized protostellar cores (e.g.. lScott & BlacMri980.) . 
Assuming that the correlation B ^ p^ exists, we predict a sim- 
ilar power-law tail in the PDF of the magnetic field strength 
with a slope from -3^ to -2| for the range of 7 e [1/2,2/3] 
and m G [-7/4,-3/2]. A power index of -2.7 measured by 
IColUns et aP (12010^ is consistent with 7 « 0.6 at p/po < 10^ 
and with their PDF slope m = -1 .64. 



5. CONCLUSIONS 

We found an intriguing agreement between the probability 
distribution of molecular cloud densities in recent observa- 
tions (Kainulainen et al. 2009) and in a deep AMR simulation 
of self-gravitating, supersonically turbulent molecular cloud 
(iPadoan et al.ll2005l) . In both cases, star-forming clouds dis- 
play strong deviations from lognormal density distribution in 
the form of power-law tails at high density. In contrast, clouds 
with no active star formation display purely lognormal distri- 
butions of density. 

We attribute the origin of the tails to the fundamental self- 
similar properties of the r~^ isothermal collapse and r"'^/^ 
pressure-free collapse laws, which control the density profiles 
of collapsing structures. This allows us to predict the power- 
law indices for the mass density {m £ [-7/4,-3/2]) and for 
the projected density {p E [-2.8,-2]) depending on the physi- 
cal conditions in the parent molecular cloud (Ms, a, etc.) that 
broadly agree with both observations and simulation results. 

Our results may suggest a reconciliation of various attempts 
to build a phenomenological theory of star formation and will 
contribute to the interpretation of numerical simulations in 
terms of the proposed p henomenologies (e.g.. ISchmidt et al.] 
lIoToHCho&KimllloTlh . 
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